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Numerical  simulations  are  carried  out  in  this  paper,  to  study  the  effect  of  air  stoichiometry  change  on  a 
single-channel  proton  exchange  membrane  (PEM)  fuel  cell  undergoing  load  change.  For  this  purpose,  the 
air  stoichiometry  change  rate  and  pattern,  as  well  as  its  lagging  time  to  the  load  change  and  initial  value,  are 
taken  into  account  in  this  model,  respectively.  The  transient  response  of  cell  potential  to  load  change  has 
been  investigated  and  the  undershoot  behavior  has  been  observed.  Detailed  results  are  further  presented 
to  show  the  transient  response  of  the  cell  in  terms  of  local  current  density  and  oxygen  mass  fraction 
distributions.  This  work  shows  that  the  extent  of  air  dilution/starvation  down  the  channel  depends  on  the 
air  stoichiomtric  flow  ratio  change  rate,  pattern,  etc.  It  is  also  shown  this  undershoot  can  be  eliminated 
by  adjusting  the  initial  air  flow  ratio,  or  by  improving  the  air  flow  ratio  change  rate,  etc.  Finally,  the 
validity  of  the  numerical  model  in  this  paper  is  verified  by  the  comparison  between  the  model  results 
and  experimental  data,  showing  that  this  numerical  model  can  be  resorted  to  reveal  the  effects  of  air 
stoichiometry  change  on  the  dynamic  behavior  of  a  PEMFC. 

©  2008  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

PEM  fuel  cell  is  regarded  as  a  potential  future  power  source  for 
electrical  vehicles,  due  to  its  high  efficiency,  zero  emission  and  fast 
start-up  at  room  temperature.  As  power  source  for  transportation, 
electrical  vehicle  often  operate  dynamically,  which  easily  cause 
some  harmful  phenomena  occurring  in  fuel  cell  [  1  -3  ].  For  example, 
during  start-up  and  acceleration  of  electrical  vehicles,  the  reac¬ 
tant  starvation  will  occur  in  fuel  cell.  In  the  case  of  fuel  starvation 
on  anode,  carbon  corrosion  and  water  electrolysis  will  occur  [2], 
according  to: 

Carbon  corrosion: 

C  +  2H20  -*  C02+4H++  46- 

Water  electrolysis: 

2H20  -*  02+4H++ 46- 

In  the  case  of  oxidant  starvation  on  cathode,  hydrogen  pump  will 
occur  [2],  according  to: 
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Hydrogen  pump: 

2H+  +  2e_  -+  H2 


The  presence  of  hydrogen  on  cathode  could  chemically  gener¬ 
ate  heat  on  the  platinum  particles  and  result  in  the  local  hot  spots 
in  the  MEA,  which  could  lead  to  failure  of  the  fuel  cells  [3].  The 
acceleration  degradation  will  shorten  the  fuel  cell  life.  So  the  main 
aim  of  transient  behavior  research  is  to  eliminate  these  harmful 
phenomena  in  fuel  cell  and  to  prolong  the  fuel  cell  life. 

Over  the  last  two  decades,  there  has  been  growing  interest 
in  studying  the  transient  behavior  of  PEMFC.  Generally,  the  pub¬ 
lished  works  in  this  field  can  be  classified  into  experimental 
studies  [4-9],  semi-empirical  modeling  [10-12]  and  mechanis¬ 
tic  modeling.  The  mechanistic  modeling  can  be  further  divided 
into  single-phase  [13-22]  and  two-phase  [23-26]  approaches.  For 
experimental  research,  Van  Zee  group  [4-6]  indicated  the  possibil¬ 
ity  of  overshoot  behavior  dependent  on  the  rate  of  voltage  change, 
and  stoichiometry.  Shen  et  al.  [9]  studied  the  effect  of  pre-supplied 
air  stoichiometry  on  cell  voltage  by  applying  a  variable  air  stoi¬ 
chiometry  during  load  change.  For  the  semi-empirical  modeling,  it 
combines  theoretically  derived  differential  and  algebraic  equations 
with  empirically  determined  relationships.  Zhang  et  al.  [10]  pro¬ 
posed  a  semi-empirical  model  based  on  the  measurements  from  a 
Nexa™  PEM  fuel  cell  module  under  different  load  conditions.  Semi- 
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Nomenclature 

a  specific  area  of  electrode  (cm-1 ) 

Ach  cross-sectional  area  of  flow  channel  (cm2) 

Amea  area  of  the  MEA  ( cm2 ) 

Ch2  local  hydrogen  concentration  (mol  m-3 ) 

Cfje2f  reference  hydrogen  concentration  (mol  m-3 ) 

Co2  local  oxygen  concentration  (mol  m-3 ) 

reference  oxygen  concentration  (mol  m-3 ) 

D  diffusion  coefficients  (m2  s-1 ) 

Ecell  cell  operating  potential  (V) 

Eoc  open-circuit  voltage  (V) 

F  Faraday’s  constant  =  96,487  (C  mol-1 ) 

ia  anode  local  current  density  (A  cm-2 ) 

ic  cathode  local  current  density  (A  cm-2 ) 

anode  reference  exchange  current  density  (A  cm-2 ) 
i£ecf  cathode  reference  exchange  current  density 

(Acm-2) 

/  cell  operating  current  density  (A  cm-2 ) 

/<p  hydraulic  permeability  ( m2 ) 

Mh2  molecular  weight  of  hydrogen  (kg  mol-1 ) 

Mh2o  molecular  weight  of  water  (kg  mol-1 ) 

Mq2  molecular  weight  of  oxygen  (kg  mol-1 ) 
nd  electro-osmotic  drag  coefficient 

p  pressure  (atm) 

R  universal  gas  constant  (8.314J  mol-1  K-1 ) 

T  temperature  (I<) 

u  velocity  vector  (ms-1) 

W  species  mass  fraction  (dimensionless) 

Greek  letters 
s  porosity 

p  overpotential 

fi  viscosity  (kg  m-1  s-1 ) 

p  density  (kg  m-3) 

a  electrode  electronic  conductivity  (S  m-1 ) 

§  stoichometric  flow  ratio 


The  objective  of  this  work  is  to  study  the  effects  of  air  stoichiom¬ 
etry  changes  on  the  transient  behavior  of  the  fuel  cell  by  applying 
an  increase  in  air  stoichiometry  during  loading,  and  the  emphasis 
is  placed  on  the  four  different  aspects  of  air  stoichiometry  change, 
including  its  rate  and  pattern,  as  well  as  its  lagging  time  to  the 
load  change  and  initial  value.  By  applying  the  2D,  isothermal,  time- 
dependent  model  proposed  in  this  work,  the  effects  of  these  four 
different  aspects  on  the  transient  characteristics  of  cell  potential 
and  mass  transport  are  investigated,  respectively. 


2.  Model  development 

2.1.  Model  assumptions 

The  assumptions  and  simplifications  adopted  in  the  present 
model  are: 

(1)  The  cell  operates  under  isothermal  conditions  at  80  °C. 

(2)  The  gas  mixture  is  an  incompressible  ideal  fluid. 

(3)  The  flow  in  the  gas  channel  is  laminar. 

(4)  The  diffusion  layer,  catalyst  layer  and  membrane  are  isotropic 
and  homogeneous,  and  the  membrane  is  impermeable  to  gas 
species. 

(5)  Ohmic  potential  drops  in  the  diffusion  layers  and  bipolar  plates 
are  neglected  due  to  their  high  electrical  conductivities. 

(6)  The  contact  resistance  between  any  two  parts  in  the  fuel  cell  is 
neglected. 

(7)  The  catalyst  layer  is  considered  as  a  thin  interface  between 
membrane  and  diffusion  layer. 

(8)  Reactant  gases  dissolving  in  the  electrolyte  at  the  catalyst  are 
not  taken  into  account. 

2.2.  Computational  domain 

Fig.  1  shows  the  schematic  diagram  of  a  PEM  fuel  cell  with  a 
single  channel  for  this  model. 


2.3.  Governing  equations 


empirical  modeling  is  very  useful  for  making  quick  predictions 
on  the  dynamics  of  fuel  cell  stacks  or  systems.  However,  it  can¬ 
not  provide  an  adequate  physical  understanding  of  the  phenomena 
inside  the  cell  [27].  For  numerical  research,  mechanistic  modeling 
has  received  much  attention  for  better  understanding  of  the  main 
transient  processes  in  PEM  fuel  cells,  including  species  transport, 
membrane  hydration/dehydration,  heat  transfer,  etc.  Among  these, 
the  species  transport  in  the  flow  channel  and  the  gas  diffusion  layer 
(GDL)  is  of  importance  on  fuel  cell  dynamics.  References  [17,18,22] 
have  studied  the  species  transport  process  and  its  effect  on  the 
dynamics  of  the  cell.  For  example,  Shimpalee  et  al.  [17]  studied  the 
effect  of  rate  change  of  cell  voltage  on  the  magnitude  of  current 
overshoot  and  found  the  peak  of  current  overshoot  depends  on  the 
rate  of  change  of  the  voltage  and  the  amount  of  gas  flow  rate.  Um 
et  al.  [18]  proposed  a  2D  transient  model  and  studied  the  effects  of 
hydrogen  dilution  along  the  anode  gas  channel.  It  is  noted  that  all 
the  predictions  except  Ref.  [22]  take  a  step  change  of  stoichiomet¬ 
ric  flow  ratio  (stoic)  or  constant  cathode/anode  flows  during  load 
change  to  study  the  transient  behavior  of  fuel  cell.  Hu  and  Fan  [22] 
take  into  account  the  effect  of  variable  flow  rate  on  the  dynamic 
behavior  of  fuel  cell  by  adopting  a  sinusoidal  change  of  air  stoi¬ 
chiometry  while  keeping  other  parameters  constant.  However,  to 
our  knowledge,  no  numerical  study  has  been  reported  by  applying 
a  variable  flow  rate  during  load  change. 


In  this  paper,  the  transient,  two-dimensional,  multi-component, 
isothermal  and  single  phase  transportation  model  can  be  written 
in  a  single  domain  forms  as  follows: 

Continuity: 

|?  +  V.(pu)  =  Sm  (1) 
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Table  1 

Source  terms  for  governing  equations  in  various  regions  of  a  PEMFC 


GDL 

CL,  a 

CL,  c 

Membrane 

Mass 

Sm  =  0 

Sm  -  Mh2Sh2  +  MH205w 

Sm  =Mh2oSo2  +Mh2o  Sw 

Sm  =  0 

Momentum 

on 

II 

1 

E4 

on 

II 

1 

nTs 

cj- 

Su  =  -jtU 

KP 

/ 

Sh2  =0 

-"ft 

1 

II 

X 

on 

o 

II 

JF 

Sh2  =0 

Species 

on 

o 

II 

o 

So2=0 

g* 

II 

1 

$l~- 

Sq2  =0 

o 

II 

Jr 

Sw  =  Nia 

CC  _  nd  ,•  1  ic 

\v  —  T lc  +  2 F 

o 

ii 

Jr 

Momentum: 


correlation  [32]. 


1 

£ 


dpu 


+  -V(puu) 


=  -Vp  +  V  ■  T  +  Su 


(2) 


/  T  \  3/2  (  P  \ 

°k  =  D°(f0)  \Pb/  f°r  §as  channel 


Species: 

^±  +  V.(uck)  =  V(DfVck)  +  Sk  (3) 


Source  terms  in  the  above  governing  equations  are  summa¬ 
rized  in  Table  1  for  various  sub-domains  of  the  fuel  cell.  In  the  gas 
channels,  the  porosity  s  becomes  unity,  and  Eq.  (2)  becomes  the 
conventional  form  of  the  momentum  equation.  In  the  porous  media 
region,  the  source  term  in  momentum  conservation  equation,  Su, 
represent  Darcy’s  drag  force. 

The  source  terms  in  Eqs.  (1  )-(3)  are  directly  related  to  the  elec¬ 
trical  current,  through  which  the  transport  equations  are  coupled 
with  electrochemical  processes.  In  this  work,  the  electrical  current 
is  modeled  by  the  Butler-Volumer  equation  [16,28,29]  as  follows: 


Anode: 

>a  =  <a( 

1/2 

exp(§F»a) 

vCH2.refy 

i 

Cathode: 

ic  =  ai'0c{  ( 

f  Co2  N 

exp 

vO>2,refy 

1 

(IN 


-  exP  N 


(§F^). 


(4) 


(5) 


where  a4,  a£,  acc  and  acA  are  the  anodic  and  cathodic  charge  transfer 
coefficient. 


2.4.  Physical  and  transport  properties 

The  water  content  and  the  flux  of  water  in  the  membrane  are 
important  for  membrane  conductivity,  which  have  a  great  influence 
on  cell  performance  and  local  current  density  distribution.  In  this 
model,  the  water  transport  equation  in  the  entire  MEA  region  is 
encompassed  in  the  general  form  of  Eq.  (3),  where  the  source/sink 
terms  are  detailed  in  Table  1 .  The  first  term  of  5^  and  5^  in  Table  1 
represents  the  electro-osmotic  drag,  and  the  second  source  term  of 
5^  is  the  water  generation  from  ORR  in  the  cathode  catalyst  layer. 
For  the  first  source  term,  nd  is  introduced  to  stand  for  the  electro- 
osmotic  drag  coefficient,  which  has  been  correlated  by  Springer  et 
al.  [30]  to  the  membrane  water  content,  A,  and  water  activity,  a. 
The  physical  parameters  of  membrane,  such  as  EW,  pdry  and  em, 
are  shown  in  Table  2. 

In  the  mass  transport  equation,  the  diffusion  coefficient  of 
species  k,  Dk,  varies  in  different  sub-regions  of  the  PEM  fuel  cell, 
depending  on  the  physical  phase  of  species  k.  In  the  anode  and  cath¬ 
ode  gas  channel,  the  diffusion  coefficient  is  calculated  as  a  function 
of  temperature  and  pressure  [31  ].  For  the  porous  regions  of  a  PEM 
fuel  cell  such  as  the  gas  diffusion  and  catalyst  layers,  the  expres¬ 
sion  is  modified  into  the  effective  species,  D£ff,  using  Bruggman 


=  sj  5Dk  for  porous  regions 

where  D0  is  the  diffusivity  in  gas  channel  at  standard  condition.  The 
mass  diffusivities  of  hydrogen  and  oxygen  through  the  membrane 
are  usually  a  few  orders  lower  than  that  in  gas  and  are  taken  to 
be  constant.  Table  2  lists  the  geometries,  physical  parameters  and 
operating  conditions  used  in  this  model. 


Table  2 

Geometries,  physical  parameters  and  operating  conditions 


Parameter 

Symbol 

Value 

Channel  length 

Ah 

20  mm 

Channel  width 

Wch 

1  mm 

Gas  diffusion  layer  thickness 

fgdl 

260  (Jim 

Membrane  thickness  (Nafion  112) 

tm 

51  [xm 

Porosity  of  anode/cathode  gas  diffusion  layers 

£gdl 

0.6 

Volume  fraction  of  ionomer  in  catalyst  layer 

£m 

0.26 

H2  reference  mole  concentration, 

CH2f 

40.88  mol  m-3 

02  reference  mole  concentration, 

40.88  mol  m-3 

Anode  exchange  current  density  multiply 
specific  area 

aC 

5.0  x  107  Am~3 

Cathode  exchange  current  density  multiply 
specific  area 

< 

1.2  x  102  Am-3 

Anode  charge  transfer  coefficient  for  anode 

0.5 

reaction 

Cathode  charge  transfer  coefficient  for  anode 

ac 

0.5 

reaction 

Anode  charge  transfer  coefficient  for  cathode 

«a 

0.5 

reaction 

Cathode  charge  transfer  coefficient  for 

ac 

0.5 

cathode  reaction 

Permeability  of  membrane 

I<m 

1.8E-18  m2 

Permeability  of  anode/cathode  gas  diffusion 

Kgdl 

2.0E-12  m2 

layers 

Dry  membrane  density 

Pdry 

2000  kg  rrr3 

Equivalent  weight  of  electrolyte  in  membrane 

EW 

1.1  kg  mol-1 

Faraday  constant 

F 

96487  C  mol”1 

Universal  constant 

R 

8.314J  mol-1  K-1 

H2  diffusivity  in  membrane 

dh2 

2.59  x  10-10  m2  s“ 

02  diffusivity  in  membrane 

dZ2 

1.22  x  10-10m2  s“ 

H2  diffusivity  in  the  anode  gas  channel  at 

D0,H2,c 

1.1  x  10-4  m2  s-1 

standard  condition 

H20  diffusivity  in  the  anode  gas  channel  at 

D0,  W,  a 

1.1  x  10-4  m2  s_1 

standard  condition, 

02  diffusivity  in  the  cathode  gas  channel  at 

D0,O2,c 

3.2  x  10-5  m2  s_1 

standard  condition 

H20  diffusivity  in  the  cathode  gas  channel  at 

Do,  w,  c 

7.3  x  10-5  m2  s”1 

standard  condition 

Cell  temperature 

T 

353 1< 

Humidity  of  anode  inlet  stream  (%) 

RH 

100 

Humidity  of  cathode  inlet  stream  (%) 

RH 

100 

Pressure  at  the  inlet  of  anode  gas  channel 

p. 

1  in,  a 

1  atm 

Pressure  at  the  inlet  of  cathode  gas  channel 

An,  c 

1  atm 

Anode  stoichiometric  flow  ratio 

2.0 
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time,  s 

Fig.  2.  Cell  average  current  density  and  cathode  stoichiometry  changes  with  time 
for  case  A. 

2.5.  Cell  potential  and  current  density 

Taking  account  the  polarization  losses,  the  actual  cell  voltage  of 
a  fuel  cell,  Ecell,  can  be  calculated  as 

£Cell  —  Eoc  ~  Pact  ~  Pohm,  mem  —  Pconc  (6) 

where  Eoc  is  the  open-circuit  voltage,  pact  the  activation  polariza¬ 
tion,  ^ohm.mem  the  ohmic  polarization  and  pConc  the  concentration 
polarization.  It  is  noted  that  the  ohmic  loss  offered  by  contact  resis¬ 
tance  is  neglected.  The  activation  losses  can  be  calculated  from  the 
Butler-Volmer  equation.  The  open-circuit  voltage  [33],  the  ohmic 
and  concentration  losses  [34]  can  be  represented  as  follows: 


y,  m 

Fig.  4.  Local  current  density  distributions  on  the  cathode  catalyst  layer  along  chan¬ 
nel  direction  for  case  Al. 

ex m  has  been  modeled  by  Springer  et  al.  [30]  in  1991.  In  the  past 
decade,  several  related  models  have  been  proposed  subsequently 
[35-37].  However,  those  models  are  quite  complicated  and  the 
physical  parameters  in  the  models  are  difficult  to  be  determined. 
In  comparison,  the  membrane  conductivity  given  by  Springer  et  al. 
is  expressed  in  water  content  and  temperature  in  a  simple  formula 
and  these  two  parameters  in  the  formula  can  be  easily  determined. 
So  this  model  is  wildly  accepted  and  used  in  most  of  the  numerical 
studies  on  PEM  fuel  cell. 

2.6.  Boundary  and  initial  conditions 


Eoc  =  0.2329  +  0.00251 

(7) 

7?ohm,mem  —  IRmcm 

(8) 

RT  /  I  \ 

(9) 

^onc  =  ^  111  (1  “  /[) 

where  T  is  the  cell  temperature,  4  is  the  limiting  current  density, 
Rmen  the  resistance  caused  by  the  electrolyte  membrane  to  the  flow 
of  H+  ions,  which  depends  on  the  thickness  of  the  membrane,  tm» 
and  the  proton  conductivity  of  the  membrane,  crm,  and  is  repre¬ 
sented  by  the  expression: 

R-mem  =  -  (10) 


2.6.1.  Inlet  boundaries 

The  inlet  values  at  the  anode  and  cathode  are  prescribed  for 
the  velocity  and  species  concentrations  (Dirichlet  boundary  condi¬ 
tions).  At  the  gas  channel  inlet,  the  inlet  velocity  uin,  is  expressed 
by  the  respective  reactant  stoichiometric  flow  ratio  of  anode  and 
cathode,  £a  and  defined  based  on  the  desired  operating  current 
density,  /,  as  follows: 


f.  I  RT  Am  1 

ln’a  ’a  2F  Pin, a  A±,a  *H,in 

(U) 

^  I  RT  Am  1 
ln’C  4F  Pin, C  7+0*0, in 

(12) 

Fig.  3.  Transient  response  of  cell  voltage  changes  with  time  for  case  A. 


y,  m 

Fig.  5.  Local  oxygen  mass  fraction  distributions  on  the  cathode  catalyst  layer  along 
channel  direction  for  case  Al. 
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Fig.  6.  Local  current  density  distributions  on  the  cathode  catalyst  layer  along  chan¬ 
nel  direction  for  case  A2. 


time,  s 

Fig.  8.  Cell  average  current  density  and  cathode  stoichiometry  changes  with  time 
for  case  B. 


where  Am  is  the  geometrical  area  of  the  membrane  and  Ach  is  the 
cross-sectional  area  of  the  gas  channel.  XH  is  molar  fraction  of 
hydrogen  and  X0  molar  fraction  of  oxygen  at  the  gas  channel  inlet. 


2.6.2.  Outlet  boundaries 

At  the  outlet  of  the  gas-flow  channels,  only  the  pressure  is  being 
prescribed  as  the  desired  electrode  pressure.  For  all  other  variables, 
fully  developed  or  no-flux  conditions  are  applied  on  the  gas  channel 
outlet  (Neumann  boundary  conditions): 


f-0,^-0 

3  n  dn 


(13) 


2.6.3.  Wall 

No-slip  and  impermeable  velocity  condition  and  no-flux  condi¬ 
tion  are  adopted 


a  =  o.*i.o,fe.o 

dn  dn 


(14) 


2.7.  Numerical  procedure 

The  time-dependent  conservation  equations  are  discretized  by 
the  finite  element  method  and  solved  by  FemLab  software.  Strin¬ 
gent  numerical  tests  are  conducted  to  achieve  the  mesh  size  and 


time  step  independence  for  the  numerical  solution.  It  is  shown 
that  the  computation  conducted  on  a  mesh  size  of  5453  elements 
with  first  time  interval  At  =  0.001  s  is  satisfactory  and  these  choices 
are  used  in  this  work  to  understand  the  transient  behavior  of  the 
reactant  gas  transport  in  PEM  fuel  cells  and  its  effect  on  cell  perfor¬ 
mance.  Due  to  the  memory  and  time  requirements  for  the  iteration 
process,  the  FemLab  software  is  run  on  Dell  Precision  workstation 
(Intel  Xeon  2.2  GHz,  1.5  GB  SDRAM). 


3.  Result  and  discussion 

In  this  section,  the  effects  of  the  cathode  stoichiometry  change 
on  the  transient  characteristics  of  a  single  cell  with  straight  channel 
to  load  change  are  studied.  In  this  study,  the  cell  average  current 
density  is  used  to  view  the  load  change  and  a  linear  increment 
from  50  to  500  mA  cm-2  in  0.1  s  is  applied  on  all  the  simulation 
cases.  Note  that  the  value  of  air  stoichiometry  in  the  whole  text  is 
given  corresponding  to  500  mA cm-2,  not  to  the  temporary  cur¬ 
rent  density.  For  the  cathode  stoichiometry,  its  change  rate  and 
pattern,  its  lagging  time  to  the  load  change  and  initial  value,  are 
investigated  and  considered  as  an  input  to  study  its  effect  on  the 
cell  performance,  respectively. 


Fig.  7.  Local  oxygen  mass  fraction  distributions  on  the  cathode  catalyst  layer  along 
channel  direction  for  case  A2. 


Fig.  9.  Transient  response  of  cell  voltage  changes  with  time  for  case  B. 
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Fig.  10.  Local  current  density  d  distributions  on  the  cathode  catalyst  layer  along 
channel  direction  for  case  B. 


Fig.  12.  Cell  average  current  density  and  cathode  stoichiometry  changes  with  time 
for  case  C. 


3.1  Effects  of  air  stoichiometry  change  rate 

Most  of  predictions  take  a  step  change  of  stoichiometric  flow 
ratio  (stoic)  or  constant  cathode/anode  flows  during  load  change 
to  study  the  transient  behavior  of  fuel  cell.  For  example,  Shimpalee 
et  al.  takes  a  constant  cathode/anode  flows  to  study  the  effects  of 
cell  voltage  change  rate  on  cell  performance,  which  corresponds  to 
an  initial  stoic  of  2.6/4.4  (2.6  for  anode  and  4.4  for  cathode  flow)  at 
0.7  V  and  a  final  stoic  of  1.2/2.4  at  0.5  V.  However,  in  practice,  i.e., 
fuel  cell  engine  for  electrical  vehicle,  during  the  vehicle  start-up  or 
acceleration,  the  hydrogen/air  flow  rates  increase  corresponding  to 
the  cell  output.  So,  to  make  the  simulations  approaching  the  prac¬ 
tical  circumstance,  it  is  necessary  to  study  the  effect  of  the  cathode 
stoichiometry  change  rate  on  the  transient  behavior  of  PEMFCs.  In 
this  study,  three  different  air  stoic  rates,  from  1.2  to  2.5  in  1  s  (case 
Al ),  in  2  s  (case  A2),  in  4  s  (case  A3),  are  inputted  to  view  its  effects 
on  cell  transient  behavior.  The  operating  current  density  and  the 
air  stoic  changes  with  the  time  are  shown  in  Fig.  2. 

Fig.  3  gives  the  transient  response  of  cell  potential  at  three 
different  air  stoic  rates.  For  the  three  potential  curves,  they  have 
the  common  starting  point  of  0.82  V  and  the  same  steady  value 
of  0.70  V,  corresponding  to  the  same  initial  and  final  conditions 
for  the  three  cases.  In  this  fig.,  it  is  clearly  seen  that  the  magni¬ 
tude  of  cell  potential  overshoot  is  affected  by  the  slope  of  air  stoic 


growth.  Lower  gas  stoic  change  rate  results  in  greater  cell  voltage 
undershoot,  lower  peak  of  cell  voltage  and  a  longer  time  to  reach 
the  steady  state.  For  case  Al,  the  cell  voltage  decreases  initially 
and  reaches  the  minimum  value  of  0.66  V  at  t=  0.28  s.  And  then  it 
increases  with  time  and  reaches  the  steady  value  at  t=  1.2  s.  When 
the  air  stoic  change  rate  is  slower  as  for  case  A2,  the  cell  voltage  has 
a  minimum  value  of  0.65  V  at  t  =  0.4  s  and  reaches  the  steady  value 
at  t=  2.2  s.  And  the  minimum  cell  voltage  is  0.64  V  at  t  =  0.5  s  and 
the  time  of  reaching  the  steady  state  is  t=4.2  s  for  case  A3. 

Figs.  4  and  5  give  the  local  current  density  and  the  local  air 
mass  fraction  profiles  on  the  cathode  catalyst  layer  along  the  chan¬ 
nel  direction  at  different  times  for  case  Al.  In  Fig.  4,  the  lower 
curve  gives  the  initial  local  current  density  distribution  with  an 
average  value  of  50  mA cm-2.  And  the  upper  four  curves  give  the 
local  current  density  distribution  after  loading,  with  an  average  of 
500  mA  cm-2.  Compared  the  upper  four  curves,  lower  and  more 
non-uniform  current  density  is  clearly  seen  downstream  along  the 
channel  direction  during  the  load  change,  particularly  under  the 
arising  of  voltage  undershoot  peak.  This  can  be  easily  explained 
by  Eq.  (5),  due  to  the  rapid  depletion  of  oxygen  from  the  inlet  to 
outlet,  shown  in  Fig.  5.  When  t>  1.1  s,  the  cell  approaches  its  final 
state,  with  more  uniform  and  higher  current  density  and  oxygen 
mass  fraction  along  the  channel  direction. 

Similar  to  the  analysis  of  case  Al,  the  simulation  cases  A2  and 
A3  are  also  analyzed.  Figs.  6  and  7  give  the  local  current  density  and 
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Fig.  11.  Local  oxygen  mass  fraction  distributions  on  the  cathode  catalyst  layer  along 
channel  direction  for  case  B. 


Fig.  13.  Transient  response  of  cell  voltage  changes  with  time  for  case  C. 
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Fig.  14.  Local  current  density  distributions  on  the  cathode  catalyst  layer  along  chan¬ 
nel  direction  for  case  C. 


the  local  air  mass  fraction  distributions  on  the  cathode  catalyst  layer 
along  the  channel  at  different  times  for  case  A2.  Compared  to  case 
Al,  case  A2  shows  more  non-uniform  current  density  distribution 
down  the  channel  during  load  change  (Fig.  6).  This  trend  is  more 
apparent  when  the  voltage  undershoot  peak  arises  (t=0.4s),  due 
to  a  lower  increasing  rate  of  air  stoichiometry  of  case  A2,  resulting 
in  a  much  heavier  reactant  dilution/starvation  down  the  channel 
(Fig.  7),  and  thus  a  substantially  larger  polarization  loss  accord¬ 
ing  to  Eq.  (5).  For  case  A3,  its  air  stoichiometry  increasing  rate  is 
the  lowest  of  the  three  cases,  resulting  in  the  heaviest  reactant 
dilution/starvation  down  the  channel  during  load  change,  and  thus 
the  largest  activation  over-potential  for  oxygen  reduction  reaction 
(ORR)  is  required  to  sustain  a  specific  operating  current  density. 

3.2.  Effects  of  lagging  time  to  load  change 

In  this  section,  the  effects  of  lagging  time  of  air  stoichiometry 
response  to  load  change  on  the  transient  behavior  is  investigated. 
For  fuel  cell  usage,  i.e.,  fuel  cell  engine  for  electrical  vehicle,  during 
start-up  or  acceleration,  the  air  compressor  need  a  response  time 
to  the  load  change.  This  time  is  called  the  lagging  time  of  the  air 
flow  to  load  change.  In  this  paper,  case  A3  without  lagging  time  is 
used  as  base  case  (case  Bl),  and  two  lagging  times,  0.5  s  (case  B2) 
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Fig.  15.  Local  oxygen  mass  fraction  distributions  on  the  cathode  catalyst  layer  along 
channel  direction  for  case  C. 
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Fig.  16.  Cell  average  current  density  and  air  stoichiometry  changes  with  time  for 
case  D. 


and  1  s  (case  B3),  are  considered  as  input  to  investigate  its  effects 
on  the  cell  performance  (seen  in  Fig.  8). 

Fig.  9  gives  the  transient  response  of  cell  potential  to  the  load 
change  at  different  level  of  lagging  time.  It  shows  that  the  peak  of 
voltage  undershoot  depends  on  the  level  of  lagging  time.  Longer 
lagging  time  results  in  greater  cell  voltage  undershoot,  lower  peak 
of  cell  voltage  and  a  longer  time  to  reach  the  final  state.  For  case  B2, 
the  cell  voltage  decreases  initially  and  reaches  the  minimum  value 
of  0.60  V  at  t=  0.74  s.  And  then  it  increases  with  time  and  reaches 
the  steady  value  at  t  =  4.7  s.  Due  to  a  longer  lagging  time  of  air  stoic 
to  load  change  for  case  B3,  the  cell  voltage  curve  shows  a  lower 
minimum  value  of  0.58  V  at  t=  1.24  s  and  reaches  the  steady  value 
at  f  =  5.2  s. 

Next,  the  local  current  density  distribution,  as  well  as  the  air 
mass  fraction  distribution  on  the  catalyst  layer  along  the  channel 
are  investigated,  as  shown  in  Figs.  10  and  11.  For  convenience,  only 
the  distributions  at  the  peak  of  cell  voltage  are  considered.  Com¬ 
parison  of  the  three  cases,  case  B3  has  the  lowest  current  density 
downstream  along  the  flow  channel,  due  to  the  much  lower  oxy¬ 
gen  mass  fraction  along  the  channel  direction.  This  is  because  case 
B3  has  a  longer  lagging  time  of  air  stoic  response  to  load  change 
than  cases  Bl  and  B2,  resulting  in  a  much  heavier  reactant  dilu¬ 
tion/starvation  down  the  channel. 


Fig.  17.  Transient  response  of  cell  voltage  changes  with  time  for  case  D. 
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Fig.  18.  Local  current  density  distributions  on  the  cathode  catalyst  layer  along  chan¬ 
nel  direction  for  case  D. 


Fig.  19.  Local  oxygen  mass  fraction  distributions  on  the  cathode  catalyst  layer  along 
channel  direction  for  case  D. 


3.3.  Effects  of  air  stoic  change  pattern  on  PEMFC  characteristics 

In  this  section,  the  effect  of  air  flow  increasing  pattern  on  the 
transient  behavior  of  PEM  fuel  cell  is  investigated.  For  example, 
for  the  fuel  cell  engine  for  transportation,  during  the  start-up  and 
acceleration  of  the  vehicle,  the  air  flow  flux  shows  a  parabolic 
growth,  corresponding  to  the  output  characteristic  of  air  compres¬ 
sor.  However,  no  numerical  or  experimental  work  on  the  effects  of 
air  stoic  increasing  form,  to  our  knowledge,  has  been  reported.  So 
in  this  paper,  a  parabolic  increasing  of  air  stoic  (case  Cl )  is  consid¬ 
ered  as  input  to  investigate  its  effect  on  the  transient  behavior  of 
fuel  cell.  For  comparison,  an  exponential  increasing  (case  C3)  is  also 
used  in  this  paper.  Fig.  12  gives  the  profiles  of  parabolic,  exponential 
and  linear  increasing  (base  case)  of  air  stoic  over  time. 

Fig.  13  gives  the  cell  voltage  variations  with  time  at  different 
shapes  of  air  stoic  growth.  It  is  seen  that  the  undershoot  phenomena 
is  almost  eliminated  for  case  Cl.  For  further  analysis,  the  local  cur¬ 
rent  density  and  air  mass  fraction  distribution  occurring  at  the  peak 
of  cell  voltage  are  also  investigated,  as  shown  in  Figs.  14  and  15.  In 
Fig.  14,  a  more  higher  and  uniform  distribution  of  local  current  den¬ 
sity  is  clearly  seen  downstream  along  the  channel  for  case  Cl,  due  to 
a  more  higher  oxygen  mass  fraction  down  the  channel  on  the  cata¬ 
lyst  layer,  shown  in  Fig.  15.  This  is  because  the  air  stoic  grows  much 
faster  for  case  Cl  than  cases  C2  and  C3,  largely  eliminating  reactant 
dilution/starvation  during  the  load  change.  Whereas,  for  case  C3, 
the  oxygen  is  almost  depleted  downstream  along  the  channel  and 
a  larger  cathode  polarization  loss  occurs. 

3.4.  Effects  of  initial  air  stoichiometry 

In  this  section,  the  effect  of  initial  air  stoichiometry  on  cell  per¬ 
formance  is  investigated.  For  this  purpose,  three  different  initial  air 
flow  ratios,  1.2  (base  case),  1.4  (case  D2)  and  1.6  (case  D3)  is  applied 
in  the  simulations  as  input  to  investigate  its  effect  on  the  transient 
behavior  of  fuel  cell.  Though  the  initial  conditions  are  different,  they 
reach  the  same  final  value  in  4  s  for  all  three  cases.  The  time  pro¬ 
files  of  average  current  density,  and  air  stoichiometry  are  shown  in 
Fig.  16. 

Fig.  17  gives  the  transient  variation  of  cell  potential  at  differ¬ 
ent  initial  air  flow  ratio.  For  the  case  D3  with  an  initial  stoic  of  1.6, 
it  has  a  peak  of  cell  potential  of  0.68  V,  approaching  to  the  steady 
value  of  the  cell.  This  shows  that  the  initial  air  stoic  has  a  great 
influence  on  the  cell  dynamic  behavior,  and  that  the  cell  voltage 
undershoot  could  be  eliminated  when  a  proper  initial  stoic  is  imple¬ 
mented.  The  local  current  density,  as  well  as  the  air  mass  fraction 


along  the  channel  is  investigated,  as  shown  in  Figs.  18  and  19.  Due 
to  the  higher  initial  value  of  air  flow  ratio,  case  D2  and  case  D3 
have  a  much  higher  reactant  concentration  on  the  catalyst  layer 
downstream  along  the  channel  direction,  largely  eliminating  the 
reactant  dilution/starvation  down  the  channel.  It  is  seen  that  the 
down-the-channel  effect  can  be  alleviated  by  applying  an  adequate 
air  stoichiometric  flow  ratio. 

4.  Model  validation 

In  this  section,  our  numerical  simulation  was  compared  with 
available  experimental  data  in  the  literature  to  validate  the  model 
proposed  in  this  paper.  As  have  been  discussed  in  the  introduction, 
many  researchers  have  studied  the  dynamic  behavior  of  PEM  fuel 
cell  experimentally  [4-9].  However,  most  of  studies  in  the  literature 
are  concentrated  on  either  step  change  of  stoichiometric  flow  ratio 
or  constant  cathode/anode  flow  rate  during  load  change  to  study 
the  transient  behavior  of  fuel  cell  [4-8].  As  far  as  we  know,  only 
Shen  et  al.  [9]  from  our  group  take  into  account  a  variable  flow  rate 
during  load  change.  So  the  experimental  data  of  Ref.  [9]  is  used  to 
verify  the  validity  of  the  proposed  model  in  this  paper. 

The  validation  process  consists  of  simulating  the  model  by 
applying  the  experimental  inputs  as  model  inputs,  and  then  com¬ 
paring  experimental  outputs  with  the  simulated  results.  In  Ref.  [9], 


Fig.  20.  Comparison  between  model  and  experimental  cell  voltage  evolution  from 
Fig.  7  in  [9]. 
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Fig.  21.  Comparison  between  model  and  experimental  cell  voltage  evolution  from 
Fig.  8  in  [9]. 

the  voltage  change  of  PEMFC  over  time  when  loading  dynamically 
from  50  to  500  mA  cm-2  in  0.5  s  was  provided  and  the  effect  of 
pre-supplied  air  stoichiometry  on  the  dynamic  behavior  was  inves¬ 
tigated.  Two  cases  were  simulated,  one  with  air  stoic  change  from 
1.3  to  2.5  in  4  s,  and  the  other  from  1.5  to  2.5  in  4  s.  The  simulated 
cell  voltage  evolutions  at  different  pre-supplied  air  stoichiometry 
are  given  in  Figs.  20  and  21  and  compared  with  the  experimental 
data.  It  is  found  that  the  model  predictions  match  the  experimen¬ 
tal  data  quite  closely.  Thus,  the  validity  of  the  numerical  model 
presented  in  this  paper  is  verified. 

5.  Conclusions 

A  2D,  isothermal,  dynamic  model  has  been  presented  to  study 
the  effect  of  air  stoichiometry  change  on  a  single-channel  poly¬ 
mer  electrolyte  fuel  cell  (PEFC).  Four  different  aspects  of  air  stoic 
change,  including  the  rate  and  pattern,  as  well  as  the  lagging  time 
to  load  change  and  the  initial  value,  are  taken  into  account  in  this 
model,  respectively.  The  cell  potential  changes  with  time  are  inves¬ 
tigated  and  the  undershoot  behavior  in  cell  potential  is  observed. 
It  is  shown  that  the  air  flow  ratio  changes  have  a  great  influence  on 
the  cell  dynamic  response  to  load  change.  Detailed  results  are  fur¬ 
ther  presented  to  show  the  transient  response  of  the  cell  in  terms 
of  local  current  density  and  oxygen  mass  fraction  distributions.  It 
is  shown  that  the  oxygen  is  progressively  depleted  downstream 
along  the  gas  channel  during  the  load  change,  resulting  in  reac¬ 
tant  gas  dilution/starvation  on  the  reaction  zone.  In  the  presence 
of  large  air  dilution,  the  over-potential  for  the  oxygen  reduction 
within  the  catalyst  layer  significantly  increases  to  sustain  a  specific 
operating  current  density,  thus  leading  to  the  cell  voltage  under¬ 
shoot  during  load  change.  It  is  also  shown  that  this  undershoot  can 
be  reduced  or  eliminated  by  improving  the  air  stoic  change  rate,  or 
by  increasing  the  initial  air  stoic  flow  ratio,  etc.  Finally,  the  validity 


of  the  proposed  model  in  this  paper  is  verified  by  the  comparison  of 
the  simulated  cell  voltage  with  the  experimental  data.  The  present 
work  is  of  importance  in  optimizing  the  control  strategy  of  PEM 
fuel  cell  to  improve  its  performance  to  load  change. 
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